Projet KIT DATA SCIENCE

Rémi Genet

Partie I - Initialisation:

Import des packages utilisés

In [1]:
import requests
import urllib.request
from bs4 import BeautifulSoup
import os
import pandas as pd
import numpy as np
from sklearn import linear_model
from sklearn import metrics
import matplotlib.pyplot as plt
import copy
from sklearn.linear_model import LassoCV
from netCDF4 import Dataset
from datetime import datetime
from datetime import timedelta
import plotly.express as px
import math
import warnings
warnings.filterwarnings('ignore')

Paramètres d'éxécution:

In [2]:
download_last_datas=False
#if selecting True the execution of some cell can be very long. In order to avoid it the results of thoose cells are
#stored in file and loaded if False is selected here
launch_long_calculations=False
#Lancement du dernier graphe dynamique, très long à executer
launch_last_graph=False
In [3]:
if launch_long_calculations:
    from global_land_mask import globe

Définition des variables globales

In [4]:
global classement_url
global glossaire_url
classement_url="https://www.vendeeglobe.org/fr/classement"
glossaire_url="https://www.vendeeglobe.org/fr/glossaire"

Partie II - Scrapping et mise en forme des données

Fonctions de scrapping et de mise en forme du dataframe

In [5]:
def get_files_names(classement_url):
    #list all the excel files of ranking in the classement page
    r=requests.get(classement_url)
    return [val.attrs.get("value") for val in BeautifulSoup(r.content.decode('utf-8'), 'html.parser').find("select",{"class":"form__input m--select onsubmit-rank"}).findAll("option") if val.attrs.get("value")!=""] if r.status_code==200 else  []

def download_all_files(files_names):
    #download all the files on the ranking pages. Delete the folder if exist and recreate it.
    delete_if_exist_and_recreate_directory(os.getcwd()+"\\downloaded_rank_files\\")
    for file in files_names:
        url='https://www.vendeeglobe.org/download-race-data/vendeeglobe_'+file+'.xlsx'
        urllib.request.urlretrieve(url, os.getcwd()+"\\downloaded_rank_files\\"+file+'.xlsx')
    return files_names

def make_dict_boat_infos(boat):
    # Function that create a dictionnaire over a boat (here a list getted with Beautifullsoup)
    boat_dict={}
    for item in boat:
        boat_dict.update({item[0][:item[0].find(":")-1]: item[0][item[0].find(":")+2:]})
        
    return boat_dict

def clean_boat_dict(list_boat_dict):
    #Reagence the dictionnaire using all the boat dictionnaire created before, create a dictionnaire 
    #where the key is the number of the boat, and each value is the other element of the dictionnaire
    #As the numéro de voile doesn't always match beetween the ranking files and the glossaire page, change the ones
    #that are different in order to match correctly
    final_dict={}
    for boat in list_boat_dict:
        nv=boat.get('Numéro de voile')
        if nv!=None and nv!='ITA 06':
            nv=nv.split(' ')[-1]
            del boat['Numéro de voile']
            try:
                nv=int(nv)
                if nv==16:
                    nv=10
            except:
                nv=int(nv[3:])
                if nv==77:
                    nv=777
        elif nv==None:
            nv=59
        else:
            nv=34
        final_dict.update({nv:boat})
    return final_dict
    
def request_boats_info(glossaire_url):
    #Function that use the different precedeing function in order to get the clean boat dictionnary
    r=requests.get(glossaire_url)
    return clean_boat_dict([make_dict_boat_infos([info.contents for info in boat.children if info!='\n']) for boat in BeautifulSoup(r.content.decode('utf-8'), 'html.parser').find("section",{"class":"boats-list"}).findAll("ul", {"class":"boats-list__popup-specs-list"})]) if r.status_code==200 else  []

def delete_if_exist_and_recreate_directory(path):
    #small function that check if a folder exist, then delete it.
    #After create it again (or for first time)
    try: 
        os.rmdir(path)
    except:
        pass
    try:
        os.mkdir(path)
    except:
        pass
    return None

def download_files():
    #Function that merge the two preceding function in order to download file more conveniently
    return download_all_files(get_files_names(classement_url))

def get_df_from_excel(file, add_extension=False):
    if add_extension:
        path=os.getcwd()+"\\downloaded_rank_files\\"+file+'.xlsx'
    else:
        path=os.getcwd()+"\\downloaded_rank_files\\"+file
    df=pd.read_excel(path, header=3)
    first_line=df[df.index==0].values
    df.columns=[first_line[0,i] if df.columns[i][:7]=="Unnamed" else df.columns[i]  for i in range(len(df.columns))]
    failer=df[df['Rang\nRank']=='RET']['Nat. / Voile\nNat. / Sail'].map(lambda x: x.split(' ')[1]).values
    df=df[df.index>0].dropna(axis=1, how="all").dropna(axis=0, how="any")
    return df,failer

def create_df_all_ranking_files():
    df, failer=get_df_from_excel(sorted(os.listdir('downloaded_rank_files'))[::-1][0])
    df["date"]=sorted(os.listdir('downloaded_rank_files'))[::-1][0][:-5]
    failers=[int(f) for f in failer]
    for ranking_file in sorted(os.listdir('downloaded_rank_files'))[::-1][1:-1]:
        df_temp, failer=get_df_from_excel(ranking_file)
        failers+=[int(f) for f in failer if int(f) not in failers]
        df_temp['date']=ranking_file[:-5]
        df=df.append(df_temp)
    return clean_df(df, failers)

def clean_df(df, failers):
    df.columns=df.columns.map(lambda x: x if x.find('\n')==-1 else x[:x.find('\n')])
    df['position']=tuple(map(convert_longitude_latitude_to_float,df['Latitude'], df['Longitude']))
    df['Latitude']=df['position'].map(lambda x: x[0])
    df['Longitude']=df['position'].map(lambda x: x[1])
    df, part_added=rename_cols(df)
    df=val_to_float(df, part_added)
    df["Nat"]=df['Nat. / Voile'].map(lambda x: x[x.find(' ')-3:x.find(' ')])
    df["Voile"]=df['Nat. / Voile'].map(lambda x: int(x[1+x.find(' '):]))
    df["Skipper"]=df["Skipper / Bateau"].map(lambda x: x[:x.find('\n')])
    df["Bateau"]=df["Skipper / Bateau"].map(lambda x: x[x.find('\n'):])
    df=df.drop(columns=['Nat. / Voile',"Skipper / Bateau"])
    df["Rang"]=df["Rang"].map(int)
    df=df[df['Voile'].map(lambda x: True if x not in failers else False)]
    df=df.set_index(df["Voile"].map(str)+df["date"])
    return df
    
def rename_cols(df):
    cols=df.columns.values
    part_added=[]
    for col in range(len(cols)):
        if "Depuis" in cols[col] :
            cols[col]="Cap "+cols[col]
    for col in range(len(cols)):
        if cols[col]=='Vitesse':
            cols[col]=cols[col]+" "+cols[col-1][4:]
            part_added.append(cols[col-1][4:])
        elif cols[col]=='VMG':
            cols[col]=cols[col]+" "+cols[col-2][4:]
        elif cols[col]=='Distance':
            cols[col]=cols[col]+" "+cols[col-3][4:]
    
    df.columns=cols
    return df,part_added

def val_to_float(df, part_added):
    #make the Vitesse, VMG, and distance as float
    for part in part_added:
        for col_type in ['Vitesse','VMG', 'Distance']:
            df[col_type+' '+part]=df[col_type+' '+part].map(lambda x: float(x[: x.find(' ')]))
    for col_type in ['DTF','DTL']:
        df[col_type]=df[col_type].map(lambda x: float(x[: x.find(' ')]))
    return df

def convert_longitude_latitude_to_float(lattitude,longitude):
    #Convert longitude and lattitude in a tupple of float
    return convert_latitude(lattitude),convert_longitude(longitude)
    
def check_if_nan_in_list(my_list):
    #function that return true if there is a Nan in a list, False if not
    #Not use in scrapping but later in the project
    return sum([1 for item in my_list if np.isnan(item)])>0

def convert_latitude(latitude):
    #formula to convert latitude in float
    d, m, s = map(float, (latitude[:2], latitude[3:5], latitude[6:8]))
    return (d + m / 60. + s / 3600.) * (1 if ('N' in latitude) else -1)

def convert_longitude(longitude):
    #formula to convert longitude in float
    d, m, s = map(float, (longitude[:2], longitude[3:5], longitude[6:8]))
    return (d + m / 60. + s / 3600.) * (1 if ('E' in longitude) else -1)

def add_boat_infos(df, boat_dict):
    #Function that create a dataframe by merging the dataframe created with ranking files and the dictionnary from the glossaire page
    df_complete=copy.deepcopy(df)
    for key in boat_dict.get(list(boat_dict.keys())[np.argmin(list(map(lambda x: len(boat_dict.get(x)), boat_dict.keys())))]).keys():
        df_complete[key]=df_complete['Voile'].map(lambda x: boat_dict.get(x).get(key))
    df_complete['Année lancement']=df_complete['Date de lancement'].map(lambda x: int(x[-4:]))
    df_complete['Longueur']=df_complete['Longueur'].map(lambda x: float(x[: x.find(' ')].replace(',','.')))
    df_complete['Largeur']=df_complete['Largeur'].map(lambda x: float(x[: x.find(' ')].replace(',','.')))
    df_complete["Tirant d'eau"]=df_complete["Tirant d'eau"].map(lambda x: float(x[: x.find(' ')].replace(',','.')))
    df_complete['Déplacement (poids)']=df_complete['Déplacement (poids)'].map(lambda x: float(x[: x.find(' ')].replace(',','.')) if x!='NC' and x!='nc' else np.nan)
    df_complete['Hauteur mât']=df_complete['Hauteur mât'].map(lambda x: float(x[: x.find(' ')].replace(',','.')))
    df_complete['Surface de voiles au près']=df_complete['Surface de voiles au près'].map(lambda x: float(x[: x.find(' ')].replace(',','.')))
    df_complete['Surface de voiles au portant']=df_complete['Surface de voiles au portant'].map(lambda x: float(x[: x.find(' ')].replace(',','.')))
    df_complete['Foil']=df_complete['Nombre de dérives'].map(lambda x: True if x in ['foils', 'foiler'] else False)
    return df_complete

Récupération des données et mise en forme dans le dataframe

In [6]:
if download_last_datas:
    download_files()
#Creation of a dataframe using the ranking file 
df=create_df_all_ranking_files()   
#Get the specific informations about the boats on the glossaire page
boat_dict=request_boats_info(glossaire_url)
#Merge all those information in a dataframe
df_complete=add_boat_infos(df, boat_dict)

Partie III - Utilisation des données au travers de modèles

A- Utilisation des données du bateau pour estimer et sa VMG

In [7]:
#We will try a regression beetween the informations we have about the boat, in our dataframe we currently have this columns
print(df_complete.columns)
Index(['Rang', 'Heure FR', 'Latitude', 'Longitude', 'Cap Depuis 30 minutes',
       'Vitesse Depuis 30 minutes', 'VMG Depuis 30 minutes',
       'Distance Depuis 30 minutes', 'Cap Depuis le dernier classement',
       'Vitesse Depuis le dernier classement',
       'VMG Depuis le dernier classement',
       'Distance Depuis le dernier classement', 'Cap Depuis 24 heures',
       'Vitesse Depuis 24 heures', 'VMG Depuis 24 heures',
       'Distance Depuis 24 heures', 'DTF', 'DTL', 'date', 'position', 'Nat',
       'Voile', 'Skipper', 'Bateau', 'Architecte', 'Chantier',
       'Date de lancement', 'Longueur', 'Largeur', 'Tirant d'eau',
       'Déplacement (poids)', 'Nombre de dérives', 'Hauteur mât',
       'Surface de voiles au près', 'Surface de voiles au portant',
       'Année lancement', 'Foil'],
      dtype='object')
In [8]:
# In order to don't have overlapping beetween speed of different classement we will 
#try to explain the ones "Depuis le dernier classement"
#We can use the informations Longueur, Largeur, Tirant d'eau, Année de lancement, Surface de voile au près et au portant
#Nombre de dérives, Hauteur Mât, and deplacement poids

#However we can see that Nombre de dérives doesn't only includes int, but also "foils" or "foiler"
print("distinct values in 'Nombre de dérives' "+str(df_complete['Nombre de dérives'].unique()))
#A solution to this is to create a dummies variables
#But as foils and foiler refer to the same things, we will use the created column Foil that contain true of false if the boat have
#foils or not
print("distinct values in 'Foil' "+str(df_complete['Foil'].unique()))
#We can also see that some informations are missing in Déplacement poids:
print("distinct values in 'Déplacement (poids)' "+str(df_complete['Déplacement (poids)'].unique()))
print("There is nan in Déplacement (poids): "+ str(check_if_nan_in_list(list(df_complete['Déplacement (poids)'].unique()))))
#A solution to this is to change the nan by the mean values of this columns
distinct values in 'Nombre de dérives' ['foils' '2' 'foiler' '2 asymétriques']
distinct values in 'Foil' [ True False]
distinct values in 'Déplacement (poids)' [8.  nan 7.6 8.5 7.8 7.7 9.  7.9 8.9 7. ]
There is nan in Déplacement (poids): True
In [9]:
#Let's first create the explained variables of our two models
Y_VMG=df_complete['VMG Depuis le dernier classement'].values.reshape(-1,1)
Y_Vitesse=df_complete['Vitesse Depuis le dernier classement'].values.reshape(-1,1)
#Let's check if both don't contain NaN values
print("There is nan in Y_VMG: "+str(check_if_nan_in_list(Y_VMG)))
print("There is nan in Y_Vitesse: "+str(check_if_nan_in_list(Y_Vitesse)))
#We can see that we don't have nan, we can go to the explanatory variables
There is nan in Y_VMG: False
There is nan in Y_Vitesse: False
In [10]:
#Transform the foil information in a dummy variable
foils=df_complete['Foil'].map(lambda x: 1 if x else 0).values
#get the weights of boats and change nan with the mean
mean_weight=np.mean([item for item in df_complete['Déplacement (poids)'].unique() if np.isnan(item)==False])
weights=df_complete['Déplacement (poids)'].map(lambda x: mean_weight if np.isnan(x) else x).values

X=np.zeros((len(weights),9 ))
X[:,-1]=foils
X[:,0]=weights
X[:,1:-1]=df_complete[['Longueur', 'Largeur', "Tirant d'eau", 'Hauteur mât',
       'Surface de voiles au près', 'Surface de voiles au portant',
       'Année lancement']]

X_colnames=['Déplacement (poids)','Longueur', 'Largeur', "Tirant d'eau", 'Hauteur mât',
       'Surface de voiles au près', 'Surface de voiles au portant',
       'Année lancement','Foil']
In [11]:
clf=linear_model.LinearRegression(fit_intercept=True, normalize=True)
clf.fit(X,Y_VMG)
print("Les coefficients de la regressions de la VMG par rapport à nos variables sont:")
print([X_colnames[i]+": "+str(clf.coef_[0,i])+", " for i in range(len(X_colnames))])

res=Y_VMG-clf.predict(X)
plt.plot(res)
plt.xlabel("observation number")
plt.ylabel("Residus")
plt.title("Resultats avec une regression linéaire classique")
plt.show()

lasso_model = LassoCV(cv=5, random_state=0,n_alphas=100, normalize=True).fit(X, Y_VMG.reshape(-1))
variable_selected_lasso=[X_colnames[var]+": "+str(lasso_model.coef_[var]) for var in range(len(lasso_model.coef_)) if np.isclose(lasso_model.coef_[var],0)==False]
residus_test_lasso=lasso_model.predict(X)-Y_VMG.reshape(-1)
print("Les coefficients de la regressions de la VMG en utilisant un modèle Lasso par rapport à nos variables sont:")
print(variable_selected_lasso)
plt.plot(residus_test_lasso)
plt.xlabel("observation number")
plt.ylabel("Residus")
plt.title("Resultats avec une regression lasso")
plt.show()
Les coefficients de la regressions de la VMG par rapport à nos variables sont:
['Déplacement (poids): -0.5148027710059779, ', 'Longueur: -0.0241040500726858, ', 'Largeur: 1.2752149920610527, ', "Tirant d'eau: -8.881784197001252e-15, ", 'Hauteur mât: -0.6762767559081085, ', 'Surface de voiles au près: 0.01714552387866285, ', 'Surface de voiles au portant: 0.0035816536034605777, ', 'Année lancement: 0.049808856797507864, ', 'Foil: 0.20314564469734672, ']
Les coefficients de la regressions de la VMG en utilisant un modèle Lasso par rapport à nos variables sont:
['Déplacement (poids): -0.5245004989948068', 'Largeur: 1.232191573232729', 'Hauteur mât: -0.6569761492547948', 'Surface de voiles au près: 0.016976307047375513', 'Surface de voiles au portant: 0.003362740678477925', 'Année lancement: 0.0479873993098772', 'Foil: 0.20207353081692594']

En regardant les résidus on peut voir que ceux-ci sont assez importants, mais surtout que ceux-ci semblent heteroscedastiques. En effet, de par la création de notre dataframe les résidus se trouvent être par date. Hors la vitesse d'un bateau ne dépend pas uniquement de ses caracteristiques mais aussi des conditions météorologiques ! Cette première regression nous permet donc d'estimer une VMG moyenne du bateau, et l'utilisation d'un modèle Lasso de trouver les variables significatives pour determiner cela. Cela ne signifie pas que les autres variables n'entrent pas en compte dans la réalité, toutefois les caracteristiques étant corrélées entre elles, prendre toutes celles-ci n'a pas d'interet.

Dans un premier temps nous pouvons donc estimer la VMG moyenne que devrait avoir chacun des bateaux en fonction de ses caracteristiques:

In [12]:
dt=df_complete['date'].max()
df_stats_boats=df_complete[df_complete['date']==dt]
foils=df_stats_boats['Foil'].map(lambda x: 1 if x else 0).values
#get the weights of boats and change nan with the mean
mean_weight=np.mean([item for item in df_stats_boats['Déplacement (poids)'].unique() if np.isnan(item)==False])
weights=df_stats_boats['Déplacement (poids)'].map(lambda x: mean_weight if np.isnan(x) else x).values

X=np.zeros((len(weights),9 ))
X[:,-1]=foils
X[:,0]=weights
X[:,1:-1]=df_stats_boats[['Longueur', 'Largeur', "Tirant d'eau", 'Hauteur mât',
       'Surface de voiles au près', 'Surface de voiles au portant',
       'Année lancement']]
VMG_boats=lasso_model.predict(X)
df_Expected_VMG=pd.DataFrame({'Voile':df_stats_boats['Voile'].values,'Rang': df_stats_boats['Rang'].values, 'VMG Expected': VMG_boats })
boats_id=df_Expected_VMG['Voile']

Nous pouvons désormais voir quel bateau réussi à battre sa VMG expected (bon choix du skipper), et ceux qui font moins bien (choix moins judicieux du skipper)

In [13]:
realized_mean_VMG=[df_complete[df_complete['Voile']==boats_id[id]]['VMG Depuis le dernier classement'].mean() for id in range(len(boats_id))]
df_Expected_realized_VMG=pd.DataFrame({'Voile':df_stats_boats['Voile'].values,'Rang': df_stats_boats['Rang'].values, 'VMG Expected': VMG_boats, 'VMG realized': realized_mean_VMG })
df_Expected_realized_VMG['Difference Realized Expected']=df_Expected_realized_VMG['VMG realized']-df_Expected_realized_VMG['VMG Expected']
In [14]:
print(df_Expected_realized_VMG.sort_values(by='Difference Realized Expected', ascending=False))
    Voile  Rang  VMG Expected  VMG realized  Difference Realized Expected
6       4     7      9.655896     12.233071                      2.577175
7      10     8      9.813134     12.145669                      2.332535
15     49    16      9.338685     11.195276                      1.856591
4      17     5     10.267844     12.122047                      1.854204
0      79     1     11.234082     12.892063                      1.657981
11      9    12     10.078128     11.653543                      1.575415
9      34    10     10.266217     11.749606                      1.483389
3      85     4     10.684474     12.152756                      1.468282
2       1     3     10.527695     11.977165                      1.449470
8    1000     9     10.410034     11.843307                      1.433273
5      18     6     11.132994     12.208661                      1.075668
1      59     2     11.234082     12.136220                      0.902138
14     99    15     10.732812     11.381890                      0.649078
16     30    17     10.484495     11.088189                      0.603694
12     53    13     10.907671     11.455118                      0.547447
10    109    11     11.321914     11.697638                      0.375724
13     27    14     11.275319     11.462992                      0.187673
26     72    27      7.657181      7.803150                      0.145969
23    777    24      8.581371      8.382677                     -0.198694
20     71    21      8.840433      8.639370                     -0.201063
21     33    22      8.771380      8.570079                     -0.201302
18     92    19     10.158195      9.934646                     -0.223549
17      7    18     11.059161     10.711024                     -0.348137
29     69    30      8.394903      7.408661                     -0.986242
19     14    20      9.783567      8.732258                     -1.051309
22      2    23     10.170506      8.653543                     -1.516963
25     83    26      9.375652      7.825984                     -1.549668
28    222    29      9.498427      7.403150                     -2.095277
27     50    28     10.166615      7.694488                     -2.472127
24     56    25     11.090971      8.077165                     -3.013806
30     11    31     10.792048      7.084252                     -3.707796
31      8    32     10.811315      6.189764                     -4.621552

En enlevant l'effet du bateau sur la VMG on peut voir que le classement ne serait alors plus le même. Cette transformation a l'avantage de permettre une comparaison plus équitable entre les skippers. Cependant nous n'en sommes encore qu'au début du Vendée Globe, hors une VMG plus faible n'est pas necessairement signe d'un moins bon skipper à l'heure actuelle, le plus important étant la VMG moyenne jusqu'à l'arrivée. En effet, il est possible et probable que certains skippers aient fait le choix d'une VMG plus faible, par exemple en s'éloignant du tracé le plus court, pour aller prendre des vents plus puissants par la suite.

Observons désormais la relation entre la difference entre réalisée et attendue entre les deux parties de la course

In [15]:
dt_sorted=np.sort(df_complete['date'].unique())
dt_first_part=dt_sorted[:int(len(dt_sorted)/2)]
dt_last_part=dt_sorted[int(len(dt_sorted)/2):]
df_complete_first_part=df_complete[df_complete['date'].map(lambda x: True if x in dt_first_part else False) ]
df_complete_last_part=df_complete[df_complete['date'].map(lambda x: True if x in dt_last_part else False) ]

realized_mean_VMG_first_part=[df_complete_first_part[df_complete_first_part['Voile']==boats_id[id]]['VMG Depuis le dernier classement'].mean() for id in range(len(boats_id))]
realized_mean_VMG_last_part=[df_complete_last_part[df_complete_last_part['Voile']==boats_id[id]]['VMG Depuis le dernier classement'].mean() for id in range(len(boats_id))]

#df_Expected_realized_VMG_last_part=pd.DataFrame({'Voile':df_stats_boats['Voile'].values,'Rang': df_stats_boats['Rang'].values, 'VMG Expected': VMG_boats last part, 'VMG realized first part': realized_mean_VMG_first_part,'VMG realized last part': realized_mean_VMG_last_part })
df_Expected_realized_VMG_by_part=pd.DataFrame({'Voile':df_stats_boats['Voile'].values,'Rang': df_stats_boats['Rang'].values, 'VMG Expected': VMG_boats, 'VMG realized first part': realized_mean_VMG_first_part, 'VMG realized last part': realized_mean_VMG_last_part })



df_Expected_realized_VMG_by_part['Difference Realized Expected first part']=df_Expected_realized_VMG_by_part['VMG realized first part']-df_Expected_realized_VMG_by_part['VMG Expected']
df_Expected_realized_VMG_by_part['Difference Realized Expected last part']=df_Expected_realized_VMG_by_part['VMG realized last part']-df_Expected_realized_VMG_by_part['VMG Expected']
df_Expected_realized_VMG_by_part['Difference Realized Expected variation']=df_Expected_realized_VMG_by_part['Difference Realized Expected last part']-df_Expected_realized_VMG_by_part['Difference Realized Expected first part']

Y=df_Expected_realized_VMG_by_part['Difference Realized Expected variation'].values.reshape(-1,1)
X=df_Expected_realized_VMG_by_part['Difference Realized Expected first part'].values.reshape(-1,1)
clf=linear_model.LinearRegression(fit_intercept=True)
clf.fit(X,Y)
df_Expected_realized_VMG_by_part
Out[15]:
Voile Rang VMG Expected VMG realized first part VMG realized last part Difference Realized Expected first part Difference Realized Expected last part Difference Realized Expected variation
0 79 1 11.234082 12.312903 13.453125 1.078821 2.219043 1.140222
1 59 2 11.234082 12.228571 12.045312 0.994489 0.811230 -0.183259
2 1 3 10.527695 12.079365 11.876562 1.551670 1.348868 -0.202803
3 85 4 10.684474 12.020635 12.282813 1.336161 1.598339 0.262178
4 17 5 10.267844 11.706349 12.531250 1.438506 2.263406 0.824901
5 18 6 11.132994 12.068254 12.346875 0.935260 1.213881 0.278621
6 4 7 9.655896 11.788889 12.670312 2.132993 3.014417 0.881424
7 10 8 9.813134 11.925397 12.362500 2.112263 2.549366 0.437103
8 1000 9 10.410034 11.576190 12.106250 1.166156 1.696216 0.530060
9 34 10 10.266217 11.184127 12.306250 0.917910 2.040033 1.122123
10 109 11 11.321914 11.736508 11.659375 0.414594 0.337461 -0.077133
11 9 12 10.078128 11.690476 11.617188 1.612348 1.539059 -0.073289
12 53 13 10.907671 11.025397 11.878125 0.117725 0.970454 0.852728
13 27 14 11.275319 10.419048 12.490625 -0.856272 1.215306 2.071577
14 99 15 10.732812 12.652381 10.131250 1.919569 -0.601562 -2.521131
15 49 16 9.338685 10.887302 11.498437 1.548617 2.159753 0.611136
16 30 17 10.484495 10.450794 11.715625 -0.033702 1.231130 1.264831
17 7 18 11.059161 10.430159 10.987500 -0.629002 -0.071661 0.557341
18 92 19 10.158195 9.574603 10.289063 -0.583592 0.130868 0.714459
19 14 20 9.783567 7.961667 9.454687 -1.821900 -0.328879 1.493021
20 71 21 8.840433 7.882540 9.384375 -0.957894 0.543942 1.501835
21 33 22 8.771380 8.285714 8.850000 -0.485666 0.078620 0.564286
22 2 23 10.170506 6.968254 10.312500 -3.202252 0.141994 3.344246
23 777 24 8.581371 7.942857 8.815625 -0.638514 0.234254 0.872768
24 56 25 11.090971 5.742857 10.375000 -5.348114 -0.715971 4.632143
25 83 26 9.375652 6.393651 9.235938 -2.982002 -0.139715 2.842287
26 72 27 7.657181 6.596825 8.990625 -1.060355 1.333444 2.393800
27 50 28 10.166615 6.582540 8.789063 -3.584075 -1.377552 2.206523
28 222 29 9.498427 6.687302 8.107812 -2.811125 -1.390614 1.420511
29 69 30 8.394903 6.388889 8.412500 -2.006014 0.017597 2.023611
30 11 31 10.792048 6.839683 7.325000 -3.952365 -3.467048 0.485317
31 8 32 10.811315 1.371429 10.932812 -9.439887 0.121497 9.561384
In [16]:
print("intercept is: "+str(clf.intercept_)+" and coefficient is: " +str(clf.coef_))
intercept is: [0.87356593] and coefficient is: [[-0.65727136]]

On peut voir sur la seconde partie de la course que l'ensemble des bateaux affichent une VMG supérieure (du fait des conditions météorologiques dans la zone où ceux-ci se trouvent). Cet effet se retrouve dans l'intercept de notre regression. Ce qu'il est interressant de voir aussi est l'effet de rattrapage, en effet, le coefficient de notre régression est négatif, cela veut dire que les bateaux ayant été le plus lentement par rapport à l'attendue sur la premiere partie sont ceux qui ont été les plus rapides sur la seconde partie de la course. On peut donc en déduire que l'hypothèse faite precedement se vérifie.

Cependant cette analyse souffre d'un biais important. Utiliser en premier lieu la VMG qui dépend du cap, et donc des choix du skipper ne donne pas un indicateur fiable de la performance des bateaux, car elle inclut fortement les choix du skipper. Afin de corriger ce biais j'ai choisi d'essayer d'estimer par bateaux la vitesse en fonction du vent et de son sens. Cela permet une comparaison plus fine des performances des bateaux en dehors des performances humaines. Ici la difficulté vient bien plus de la récupération de la donnée, et du raccrochement à nos données initiales. Le code ci-dessous se charge de recupérer les données et de les ajouter à notre dataframe. Les données proviennent du NOAA (National Oceanic and Atmospheric Administration), de deux fichiers en format .nc dont l'un contient la vitesse du vent méridionnale et le second du vent zonal. Ces données sont observées 4x par jour sur une grille de longitude et de lattitude. De ces deux valeurs en m.s-1 est recalculé un angle du vent dans le même référentiel que celui utilisé pour le cap du bateau, ainsi que sa force. 3 colonnes sont alors ajoutées contenant l'angle du vent, l'angle du vent par rapport au bateau et sa vitesse, et ce à chaque date en fonction des données les plus proches géographiquement et temporellement.

B- Utilisation des données de vents pour estimer la VMG

In [17]:
def download_wind_data_files():    
    url_current_year_Vvalues="ftp://ftp2.psl.noaa.gov/Datasets/ncep.reanalysis/surface_gauss/vwnd.10m.gauss.2020.nc"
    url_last_year_Vvalues="ftp://ftp2.psl.noaa.gov/Datasets/ncep.reanalysis/surface_gauss/vwnd.10m.gauss.2019.nc"
    url_current_year_Uvalues="ftp://ftp2.psl.noaa.gov/Datasets/ncep.reanalysis/surface_gauss/uwnd.10m.gauss.2020.nc"
    url_last_year_Uvalues="ftp://ftp2.psl.noaa.gov/Datasets/ncep.reanalysis/surface_gauss/uwnd.10m.gauss.2019.nc"
    delete_if_exist_and_recreate_directory(os.getcwd()+"\\wind_datas\\")
    urllib.request.urlretrieve(url_current_year_Vvalues, os.getcwd()+"\\wind_datas\\current_year_Vvalues.nc")
    urllib.request.urlretrieve(url_last_year_Vvalues, os.getcwd()+"\\wind_datas\\last_year_Vvalues.nc")
    urllib.request.urlretrieve(url_current_year_Uvalues, os.getcwd()+"\\wind_datas\\current_year_Uvalues.nc")
    urllib.request.urlretrieve(url_last_year_Uvalues, os.getcwd()+"\\wind_datas\\last_year_Uvalues.nc")
    
def get_lon_lat_time_from_nc(nc_file):
    #function returning the list of available longitude and latitude in the file
    variables_in_files=nc_file.variables
    return variables_in_files['lon'][:].data,variables_in_files['lat'][:].data,variables_in_files["time"][:].data
        
def get_values_wind(Unc_file, Vnc_file, lat_in_file, lon_in_file, time_in_file, real_lat, real_lon, real_time ):
    #Using a real longitude, latitude and time wanted, return the values the closest in the space and time from the nc file
    lat_index=get_closest_index(lat_in_file, real_lat)
    if real_lon>0:
        real_lon=real_lon
    else:
        real_lon=360+real_lon
    lon_index=get_closest_index(lon_in_file, real_lon)
    time_index=get_closest_index(time_in_file, real_time)
    return calculate_wind_speed_angle(Unc_file.variables['uwnd'][time_index, lat_index,lon_index].data, Vnc_file.variables['vwnd'][time_index,lat_index,lon_index].data)

def calculate_wind_speed_angle(Uwind, Vwind):
    #Use zonal and meridional wind to calculate a wind speed and angle
    wind_speed=np.sqrt(Uwind**2+Vwind**2)
    if wind_speed!=0:
        wind_angle=put_360(math.degrees(np.arccos(-Vwind/wind_speed)))
    else:
        wind_angle=0
    return wind_speed, wind_angle

def put_360(degrees_180):
    #convertion function
    if degrees_180<0:
        degrees_180=360-degrees_180
    return degrees_180

def get_closest_index(available_values, searched_value):
    #return the index of the closest value in a list of a searched value
    return np.argmin(abs(available_values-searched_value))

def transform_date(date_df):
    y=int(date_df[:4])
    m=int(date_df[4:6])
    d=int(date_df[6:8])
    h=int(date_df[9:11])
    return (datetime(year=y, month=m, day=d, hour=h)-datetime(2020,1,1)).total_seconds()/3600

if download_last_datas:
    download_wind_data_files()

cy_Vwind = Dataset(os.getcwd()+"\\wind_datas\\current_year_Vvalues.nc", "r", format="NETCDF4")
cy_Uwind = Dataset(os.getcwd()+"\\wind_datas\\current_year_Uvalues.nc", "r", format="NETCDF4")
ly_Vwind = Dataset(os.getcwd()+"\\wind_datas\\last_year_Vvalues.nc", "r", format="NETCDF4")
ly_Uwind = Dataset(os.getcwd()+"\\wind_datas\\last_year_Uvalues.nc", "r", format="NETCDF4")
lon_in_file,lat_in_file, time_in_file_cy=get_lon_lat_time_from_nc(cy_Vwind)
time_in_file_cy-=time_in_file_cy[0]
_, _, time_in_file_ly=get_lon_lat_time_from_nc(ly_Uwind)
time_in_file_ly-=time_in_file_ly[0]

df_complete['time_delta_start_year']=df_complete['date'].map(transform_date)
real_lats=df_complete['Latitude'].values
real_lons=df_complete['Longitude'].values
real_times=df_complete['time_delta_start_year'].values
wind_speed=[]
wind_angle=[]
for iter in range(len(real_times)):
    vspeed, vangle=get_values_wind(cy_Uwind,cy_Vwind, lat_in_file, lon_in_file, time_in_file_cy, real_lats[iter],real_lons[iter],real_times[iter])
    wind_speed.append(vspeed)
    wind_angle.append(vangle)
    
    
    
df_complete['Wind speed']=wind_speed
df_complete['wind angle']=wind_angle
df_complete['wind angle with boat']=df_complete['wind angle']-df_complete['Cap Depuis le dernier classement'].map(lambda x: int(x[:-1]))
df_complete['wind angle with boat']=df_complete['wind angle with boat'].map(lambda x: abs(x) if abs(x)<180 else abs(x)-90).map(lambda x: x if x<180 else abs(x-180)).map(lambda x: abs(x-180))

print(df_complete[df_complete['Rang']==1][df_complete['date']=="20201116_170000"]['Wind speed'])
print(df_complete[df_complete['Rang']==1][df_complete['date']=="20201116_170000"]['wind angle'])
9920201116_170000    7.56439
Name: Wind speed, dtype: float64
9920201116_170000    65.806789
Name: wind angle, dtype: float64
In [18]:
df_complete
Out[18]:
Rang Heure FR Latitude Longitude Cap Depuis 30 minutes Vitesse Depuis 30 minutes VMG Depuis 30 minutes Distance Depuis 30 minutes Cap Depuis le dernier classement Vitesse Depuis le dernier classement ... Nombre de dérives Hauteur mât Surface de voiles au près Surface de voiles au portant Année lancement Foil time_delta_start_year Wind speed wind angle wind angle with boat
7920201129_110000 1 11:30 FR\n -38.851667 7.113889 66° 22.1 15.9 11.0 61° 19.1 ... foils 29.0 350.0 560.0 2019 True 8003.0 1.603122 3.576317 122.576317
5920201129_110000 2 10:30 FR\n-60min -38.258889 0.436667 51° 19.3 9.4 9.6 51° 18.4 ... foils 29.0 350.0 560.0 2019 True 8003.0 5.147815 97.815296 133.184704
120201129_110000 3 11:30 FR\n -37.900278 -0.467500 126° 16.5 16.1 8.3 136° 15.8 ... 2 28.0 300.0 620.0 2007 False 8003.0 3.360060 126.528858 170.528858
8520201129_110000 4 11:30 FR\n -38.475000 -0.991111 49° 15.6 7.1 7.8 77° 13.8 ... foils 27.4 300.0 600.0 2010 True 8003.0 5.220153 106.699241 150.300759
1720201129_110000 5 11:30 FR\n -37.386111 -1.171389 55° 14.8 7.8 7.4 52° 18.0 ... foils 29.0 310.0 550.0 2015 True 8003.0 3.360060 126.528858 105.471142
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
7220201108_140000 27 15:30 FR\n 46.439722 -1.825556 236° 10.9 10.8 0.2 357° 0.0 ... 2 29.0 260.0 580.0 1998 False 7502.0 5.923681 158.198588 71.198588
2720201108_140000 28 15:26 FR\n-4min 46.427222 -1.806111 238° 13.8 13.5 0.2 358° 0.0 ... foils 27.0 300.0 650.0 2007 True 7502.0 5.923681 158.198588 70.198588
420201108_140000 29 15:29 FR\n-1min 46.437500 -1.820278 235° 13.4 13.4 0.7 357° 0.0 ... foiler 29.0 260.0 600.0 2019 True 7502.0 5.923681 158.198588 71.198588
5020201108_140000 30 15:28 FR\n-2min 46.427500 -1.809444 237° 11.4 11.3 0.4 358° 0.0 ... 2 28.0 330.0 600.0 2006 False 7502.0 5.923681 158.198588 70.198588
22220201108_140000 31 15:30 FR\n 46.434722 -1.805833 234° 12.1 12.0 0.2 358° 0.0 ... 2 28.0 270.0 580.0 2007 False 7502.0 5.923681 158.198588 70.198588

4060 rows × 41 columns

La première approche, intuitive serait de regresser uniquement la vitesse et l'angle du vent sur le bateau. Cependant les résultats sont évidement assez mauvais, car la force appliquée par le vent n'est pas linéaire de ces deux variables. Un dataframe temporaire est ici créé avec des variables tel que le sinus et le cosinus de l'angle du vent. Enfin le bateau n'avancant pas sans vent, contrairement à la plupart des modèles il faut ici bien enlever l'intercept qui n'aurait alors pas de sens. Enfin les données sont réduites en séparant le dataframe par bateau (100 points à l'heure actuelle), et les skippers ne se mettant pas face au vent certaines configurations ne sont pas représentées ou dans des zones de vent nul. Pour contrer cela, j'ai agrandi les données d'étude en y ajoutant des valeurs dans ces cas-ci, afin de forcer le modèle à respecter ces contraintes. Enfin pour utiliser les interactions entre les variables et les variables au carré j'ai fait appel à PolynomialFeatures de sklearn.

In [19]:
from sklearn.preprocessing import PolynomialFeatures
In [20]:
boats_numbers=df_complete['Voile'].unique()
df_for_reg=copy.deepcopy(df_complete)
dt_part=dt_sorted[3:-3]
df_for_reg=df_for_reg[df_for_reg['date'].map(lambda x: True if x in dt_part else False)]
df_for_reg['div']=df_for_reg['wind angle with boat'].map(lambda x: 1/np.log(x+2))

df_for_reg['sin']=df_for_reg['wind angle with boat'].map(lambda x: np.sin(math.radians(x)))
df_for_reg['cos']=df_for_reg['wind angle with boat'].map(lambda x: np.cos(math.radians(x)))

clf=linear_model.LinearRegression(fit_intercept=False)
R2=[]
boats_model_lasso={}
boats_model_clf={}
R22=[]
test=[]
poly = PolynomialFeatures(degree=2,include_bias=False, interaction_only=False)

for boat in boats_numbers:
    Y=df_for_reg[df_for_reg['Voile']==boat][df_for_reg['Vitesse Depuis le dernier classement']>0.01]['Vitesse Depuis le dernier classement']
    X=df_for_reg[df_for_reg['Voile']==boat][df_for_reg['Vitesse Depuis le dernier classement']>0.01][['wind angle with boat', 'Wind speed', 'div','cos','sin']]
    Xbis=np.zeros((int(np.shape(X)[0]*1.3),np.shape(X)[1]))
    Xbis[:np.shape(X)[0],:]=X
    Xbis[np.shape(X)[0]:,1]=0
    Xbis[np.shape(X)[0]:,0]=np.random.randint(0,180,np.shape(Xbis)[0]-np.shape(X)[0])
    Xbis[np.shape(X)[0]:,2]=1/np.log(np.random.randint(0,180,np.shape(Xbis)[0]-np.shape(X)[0])+2)
    Xbis[np.shape(X)[0]:,3]=np.sin(np.radians(np.random.randint(0,180,np.shape(Xbis)[0]-np.shape(X)[0])))
    Xbis[np.shape(X)[0]:,4]=np.cos(np.radians(np.random.randint(0,180,np.shape(Xbis)[0]-np.shape(X)[0])))
    Ybis=np.zeros((len(Xbis)))
    Ybis[:np.shape(X)[0]]=Y
    Ybis[np.shape(X)[0]:]=0
    #Xbis=X.values
    #Ybis=Y.values
    Xbis=poly.fit_transform(Xbis)
    lasso_model=LassoCV(fit_intercept=False, normalize=True).fit(Xbis, Ybis)
    #lasso_model = linear_model.Lasso(max_iter=100000, selection='random',fit_intercept=False, normalize=False).fit(Xbis, Ybis)
    clf=linear_model.LinearRegression(fit_intercept=False).fit(Xbis,Ybis)
    pred2=lasso_model.predict(Xbis)
    pred=clf.predict(Xbis)
    #res=Y-pred
    #res2=Y-pred2
    plt.plot(range(len(pred)),pred, label='predit')
    plt.plot(range(len(pred)),pred2, label='predit Lasso')
    plt.plot(range(len(pred)),Ybis, label='True')
    R2.append(metrics.r2_score(Ybis, pred))
    R22.append(metrics.r2_score(Ybis, pred2))
    plt.legend()
    plt.title("Valeurs réelle vs estimé pour le bateau numéro "+str(boat))
    plt.show()
    boats_model_clf.update({boat: clf})
    boats_model_lasso.update({boat: lasso_model})
    test.append((np.min(X.values[:,1]),np.max(X.values[:,1]), np.mean(X.values[:,1])))

print("R2 moyen obtenu: "+str(np.mean(R2)))
print("R2 moyen obtenu: "+str(np.mean(R22)))
R2 moyen obtenu: 0.7665506250489773
R2 moyen obtenu: 0.16009505525106243
In [ ]:
 

On peut voir que les résultats de la regression, sans être totalement parfaits semblent satisfaisants. Il était dur d'esperer beaucoup plus étant donné d'une part l'approximation des données au niveau géographique et temporel, mais aussi d'autres données météorologiques pouvant ralentir la navigation comme la houle, les précipitations, ou les changements de cap opérés par le skipper. Cependant ces résultats me semblent suffisants pour pouvoir comparer les bateaux entre eux. Une première question qui nous était proposée était d'estimer l'effet des foils sur les bateaux. Je ne l'avais pas fait pour la VMG, car cela n'aurait pas eu beaucoup de sens, mais ici cela devient bien plus interressant. Je vais donc entamer la comparaison entre les bateaux pour des conditions de vents différentes. On voit ici qu'il existe tout de même une erreur importante notamement sur les zones de vent nul, avec des erreurs assez importantes.

Pour ce faire il est necessaire de prendre en compte les autres effets sur le bateaux en plus des foils. Nous allons donc regresser la vitesse prédite pour chaque bateau et dans des conditions de vents differentes sur leurs caractéristiques.

In [21]:
def transform_data_for_model(ws, wa):
    watransform=np.log(wa+2)
    return np.array([wa, ws,1/watransform, np.sin(math.radians(wa)), np.cos(math.radians(ws))]).reshape(1,-1)




df_infos_per_boat_only=df_for_reg[df_for_reg['date'].map(lambda x: True if x in dt_part[-1] else False)]
df_infos_per_boat_only=df_infos_per_boat_only[['Voile','Longueur', 'Largeur', "Tirant d'eau", 'Hauteur mât',
       'Surface de voiles au près', 'Surface de voiles au portant',
       'Année lancement','Foil']]
nb_boat_with_foil=df_infos_per_boat_only['Foil'].sum()
nb_boat_without_foil=len(df_infos_per_boat_only)-df_infos_per_boat_only['Foil'].sum()
wind_speed_step=0.25
wind_speed_max=18
wind_speed_min=0
wind_angles=[0,45,90,135,180]
speed_boat_with_foil=np.zeros((int((wind_speed_max-wind_speed_min)/wind_speed_step), len(wind_angles)))
speed_boat_without_foil=np.zeros((int((wind_speed_max-wind_speed_min)/wind_speed_step), len(wind_angles)))
Y=np.zeros((int((wind_speed_max-wind_speed_min)/wind_speed_step*len(boats_numbers))))
X=np.zeros((int((wind_speed_max-wind_speed_min)/wind_speed_step*len(boats_numbers)), 8))
windangle=90
col_names=['Longueur', 'Largeur', "Tirant d'eau", 'Hauteur mât',
       'Surface de voiles au près', 'Surface de voiles au portant',
       'Année lancement','Foil']
col=0
for windangle in wind_angles:
    count=0
    for boat in boats_numbers:
        clf=boats_model_lasso.get(boat)
        line=0
        for windspeed in np.arange(wind_speed_min, wind_speed_max, wind_speed_step):
            Y[count]=clf.predict(poly.fit_transform([transform_data_for_model(windspeed, windangle)[0]]))
            X[count,:]=df_infos_per_boat_only[df_infos_per_boat_only['Voile']==boat].values[0,1:]
            if X[count,-1]:
                speed_boat_with_foil[line, col]+=Y[count]/nb_boat_with_foil
            else:
                speed_boat_without_foil[line, col]+=Y[count]/nb_boat_without_foil
            line+=1
            count+=1
    lasso_model = LassoCV(cv=5, random_state=0,n_alphas=100, fit_intercept=False, normalize=True).fit(X, Y)
    print("les variables du modeles lasso pour un angle de "+str(windangle)+" degré sont :")
    print(lasso_model.coef_)
    print("les variables selectionnés sont ainsi:")
    print([col_names[var] for var in range(len(lasso_model.coef_)) if np.isclose(lasso_model.coef_[var],0)==False])
    col+=1

for col in range(len(wind_angles)):
    plt.plot(np.arange(wind_speed_min, wind_speed_max, wind_speed_step),speed_boat_with_foil[:,col], label='bateaux avec foil')
    plt.plot(np.arange(wind_speed_min, wind_speed_max, wind_speed_step),speed_boat_without_foil[:,col], label='bateaux sans foil')
    plt.legend()
    plt.xlabel("Vitesse du vent (m/s)")
    plt.ylabel("Vitesse du navire (noeuds)")
    plt.title("resultat avec un angle d'incidence du vent de "+str(wind_angles[col])+" degrés")
    plt.show()
les variables du modeles lasso pour un angle de 0 degré sont :
[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 6.25104997e-21 0.00000000e+00]
les variables selectionnés sont ainsi:
[]
les variables du modeles lasso pour un angle de 45 degré sont :
[ 0.          0.          0.          0.          0.         -0.
  0.00397533  0.        ]
les variables selectionnés sont ainsi:
['Année lancement']
les variables du modeles lasso pour un angle de 90 degré sont :
[ 0.          0.          0.          0.          0.         -0.
  0.00796056  0.        ]
les variables selectionnés sont ainsi:
['Année lancement']
les variables du modeles lasso pour un angle de 135 degré sont :
[0.         0.         0.         0.         0.         0.
 0.01180511 0.        ]
les variables selectionnés sont ainsi:
['Année lancement']
les variables du modeles lasso pour un angle de 180 degré sont :
[0.         0.         0.         0.         0.         0.
 0.01539506 0.        ]
les variables selectionnés sont ainsi:
['Année lancement']

Bien que le modèle soit loin de présenter une précision absolue dans le sens où de nombreuses conditions sont ignorées, on peut voir que les résultats montrent tout de même une certaine cohérence sur certains points. Les courbes partent toutes d'environ 0, en vent nul le bateau n'avance donc presque pas. Quand celui-ci augmente sa vitesse augmente et ce en fonction de l'angle du vent. Face au vent de même le bateau n'avance presque pas. Enfin il est interressant que bien qu'un lasso ne détecte pas l'effet des foils directement mais l'année de lancement du bateau comme explicatif de la vitesse, l'analyse graphique des courbes de vitesses prédites permet de voir que les bateaux avec foils sont légèrement avantagés sur les autres. On peut d'ailleurs voir que les bateaux possédant des foils sont principalement les bateaux récents. L'information année de lancement contient donc en partie l'information foil, mais aussi d'autres informations expliquant sa séléction à contrario de la variable foil.

Lien entre Vitesse et VMG

In [22]:
Y=df_complete['VMG Depuis le dernier classement'].values.reshape(-1,1)
X=df_complete["Vitesse Depuis le dernier classement"].values.reshape(-1,1)
clf=linear_model.LinearRegression(fit_intercept=False).fit(X,Y)
pred=clf.predict(X)
plt.plot(Y, label='True')
plt.plot(pred, label='predit')
plt.legend()
plt.show()
print("Intercept: "+str(clf.intercept_))
print("coefficient "+str(clf.coef_))
print("R2 de la regression: "+str(metrics.r2_score(pred,Y)))
Intercept: 0.0
coefficient [[0.85829144]]
R2 de la regression: 0.5864061334270863

Le lien entre vitesse et VMG est assez évident mais il apparait tout de même utile de le verifier, on voit ici que le modèle fonctionne très bien avec un R2 important. Ces deux variables auraient d'ailleurs pu être moins liées dans le cas où les routes prises par les skippers auraient fortement différé, auquel cas une vitesse plus rapide aurait pu entrainer une VMG plus petite car fesant un détour important. Ce n'est ici pas le cas, cela nous permet donc de voir que la vitesse des bateaux est ici un élément déterminant dans la VMG. La VMG impliquant méchaniquement le classement, l'observation de la vitesse théorique de chaque bateau a donc bien eu un interêt certain pour comprendre le classement.

Ayant désormais le lien entre la VMG et la vitesse, et ayant un modèle, bien qu'imparfait, de calcul théorique de la vitesse pour chaque bateau il est possible de recréer la VMG attendue pour chaque bateau dans des conditions similaires à celles rencontrées sur la course. Pour ce faire j'ai récupéré les données de vents de tous les bateaux durant la course. La fonction get_expected_vmg_boat calcule alors la VMG du bateau pour ces données en fonction de la vitesse prédite de celui-ci, et en fait la moyenne pour l'ensemble des données de vent.

In [23]:
wind_speed_datas=df_complete['Wind speed'].values
wind_angle_datas=df_complete['wind angle with boat'].values
In [24]:
def get_expected_vmg_boat(clf_vmg_vitesse, boat_speed_model, winds_speed_datas, wind_angle_datas):
    vmg_e=[]
    for iter in range(len(wind_speed_datas)):
        vmg_e.append(clf_vmg_vitesse.predict(boat_speed_model.predict(poly.fit_transform([transform_data_for_model(winds_speed_datas[iter], wind_angle_datas[iter])[0]])).reshape(-1,1))[0][0])
    return np.mean(vmg_e)
In [25]:
Expected_VMG=[]
for boat in boats_numbers:
    #Expected_VMG.append(clf.predict(boats_model.get(boat).predict(poly.fit_transform([transform_data_for_model(mean_wind, mean_wind_angle)[0]])).reshape(-1,1))[0][0])
    Expected_VMG.append(get_expected_vmg_boat(clf, boats_model_lasso.get(boat), wind_speed_datas,wind_angle_datas))
df_Expected_VMG=pd.DataFrame({'Voile':df_stats_boats['Voile'].values,'Rang': df_stats_boats['Rang'].values, 'VMG Expected': VMG_boats, 'New VMG expected': Expected_VMG })


realized_mean_VMG=[df_complete[df_complete['Voile']==boats_id[id]]['VMG Depuis le dernier classement'].mean() for id in range(len(boats_id))]
Foil=[df_complete[df_complete['Voile']==boats_id[id]]['Foil'].values[0] for id in range(len(boats_id))]
Annee=[df_complete[df_complete['Voile']==boats_id[id]]['Année lancement'].values[0] for id in range(len(boats_id))]

df_Expected_realized_VMG=pd.DataFrame({'Voile':df_stats_boats['Voile'].values,'Rang': df_stats_boats['Rang'].values, 'VMG Expected': VMG_boats, 'New VMG expected': Expected_VMG, 'VMG realized': realized_mean_VMG, 'Foil': Foil, 'Année lancement':Annee })
df_Expected_realized_VMG['Difference Realized Expected']=df_Expected_realized_VMG['VMG realized']-df_Expected_realized_VMG['VMG Expected']
df_Expected_realized_VMG['Difference Realized New Expected']=df_Expected_realized_VMG['VMG realized']-df_Expected_realized_VMG['New VMG expected']
df_Expected_realized_VMG['Difference Expected New Expected']=df_Expected_realized_VMG['VMG Expected']-df_Expected_realized_VMG['New VMG expected']
In [26]:
df_Expected_realized_VMG
Out[26]:
Voile Rang VMG Expected New VMG expected VMG realized Foil Année lancement Difference Realized Expected Difference Realized New Expected Difference Expected New Expected
0 79 1 11.234082 10.979083 12.892063 True 2019 1.657981 1.912980 0.254999
1 59 2 11.234082 8.044460 12.136220 True 2019 0.902138 4.091760 3.189622
2 1 3 10.527695 8.548866 11.977165 False 2007 1.449470 3.428299 1.978829
3 85 4 10.684474 10.702628 12.152756 True 2010 1.468282 1.450128 -0.018154
4 17 5 10.267844 9.933646 12.122047 True 2015 1.854204 2.188401 0.334197
5 18 6 11.132994 9.112709 12.208661 True 2015 1.075668 3.095952 2.020285
6 4 7 9.655896 10.447120 12.233071 True 2019 2.577175 1.785951 -0.791224
7 10 8 9.813134 10.253091 12.145669 True 2015 2.332535 1.892578 -0.439957
8 1000 9 10.410034 8.753782 11.843307 False 2008 1.433273 3.089525 1.656252
9 34 10 10.266217 8.654377 11.749606 True 2015 1.483389 3.095230 1.611841
10 109 11 11.321914 9.285621 11.697638 True 2010 0.375724 2.412017 2.036292
11 9 12 10.078128 8.884127 11.653543 False 2007 1.575415 2.769417 1.194002
12 53 13 10.907671 8.228561 11.455118 False 2007 0.547447 3.226558 2.679111
13 27 14 11.275319 8.100679 11.462992 True 2007 0.187673 3.362313 3.174640
14 99 15 10.732812 9.899668 11.381890 True 2019 0.649078 1.482222 0.833144
15 49 16 9.338685 8.384313 11.195276 False 2007 1.856591 2.810963 0.954372
16 30 17 10.484495 8.256769 11.088189 False 2011 0.603694 2.831420 2.227726
17 7 18 11.059161 8.162748 10.711024 True 2007 -0.348137 2.548276 2.896413
18 92 19 10.158195 9.941459 9.934646 True 2007 -0.223549 -0.006813 0.216736
19 14 20 9.783567 8.318012 8.732258 True 2007 -1.051309 0.414246 1.465555
20 71 21 8.840433 8.364689 8.639370 False 2007 -0.201063 0.274681 0.475744
21 33 22 8.771380 7.790242 8.570079 False 2000 -0.201302 0.779836 0.981138
22 2 23 10.170506 7.430668 8.653543 True 2020 -1.516963 1.222875 2.739838
23 777 24 8.581371 6.413031 8.382677 False 1999 -0.198694 1.969646 2.168340
24 56 25 11.090971 7.858289 8.077165 True 2015 -3.013806 0.218877 3.232683
25 83 26 9.375652 5.751006 7.825984 False 2006 -1.549668 2.074978 3.624646
26 72 27 7.657181 6.261177 7.803150 False 1998 0.145969 1.541972 1.396003
27 50 28 10.166615 6.278842 7.694488 False 2006 -2.472127 1.415646 3.887773
28 222 29 9.498427 6.277156 7.403150 False 2007 -2.095277 1.125994 3.221271
29 69 30 8.394903 5.253580 7.408661 False 2005 -0.986242 2.155081 3.141323
30 11 31 10.792048 7.007792 7.084252 True 2019 -3.707796 0.076459 3.784255
31 8 32 10.811315 8.659584 6.189764 True 2018 -4.621552 -2.469821 2.151731

En recréant le premier tableau on voit ici que la différence entre le réalisé présentent un biais haussier, dû en parti au fait que le modèle n'inclut pas d'intercept ici contrairement au premier. Le modele n'est donc pas parfait, durant mes recherches d'ailleurs en utilisant d'autres combinaisons de variables, et en utilisant des intercepts les résultats ici semblaient plus interressant dans le sens où ils étaient en absolu plus proche des valeurs réalisé et moins biaiser. Cependant ces modèle présentaient alors moins de sens car présentant des vitesses pouvant être positive en vent nul, et pouvant passer négative en cas de vent trop important.

In [27]:
print("sur l'ensemble des bateaux")
print("moyenne realisé: "+str(df_Expected_realized_VMG['VMG realized'].mean()))
print("moyenne expected: "+str(df_Expected_realized_VMG['VMG Expected'].mean()))
print("moyenne new expected: "+str(df_Expected_realized_VMG['New VMG expected'].mean()))
sur l'ensemble des bateaux
moyenne realisé: 10.140794497462009
moyenne expected: 10.141162594332261
moyenne new expected: 8.31993054771115

En regardant sur l'ensemble des bateaux, evidemment la premiere facon de calculer la VMG est par construction égale à celle réalisé (regression avec intercept), et différentes pour celle que nous vennons de créer.

In [28]:
print("sur les 10 premiers")
print("moyenne realisé: "+str(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']<10]['VMG realized'].mean()))
print("moyenne expected: "+str(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']<10]['VMG Expected'].mean()))
print("moyenne new expected: "+str(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']<10]['New VMG expected'].mean()))
sur les 10 premiers
moyenne realisé: 12.190106792206528
moyenne expected: 10.551137164637618
moyenne new expected: 9.641709561856329

On voit que l'en ne regardant plus que les 10 premiers bateaux la nouvelle formule permet de calculer une VMG semblent assez similaire aux valeurs obtenu à l'aide de l'ancienne formule.

In [29]:
print("sur les 10 derniers (hors bateau 8)")
print("moyenne realisé: "+str(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']>21][df_Expected_realized_VMG['Rang']<32]['VMG realized'].mean()))
print("moyenne expected: "+str(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']>21][df_Expected_realized_VMG['Rang']<32]['VMG Expected'].mean()))
print("moyenne new expected: "+str(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']>21][df_Expected_realized_VMG['Rang']<32]['New VMG expected'].mean()))
sur les 10 derniers (hors bateau 8)
moyenne realisé: 7.8903149606299205
moyenne expected: 9.449905462264422
moyenne new expected: 6.6321784475264165

Sur la fin de peloton toutefois les valeurs sont bien plus faible qu'avec l'ancienne formule. L'avantage du modele pour comparer les competences est ici clair, il n'y a plus de biais en fonction de la VMG réalisé appliqué par le modèle qui pénalise mécaniquement les derniers.

In [30]:
plt.plot(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']<32]['Difference Realized Expected'], label='ecart à expected')
plt.plot(df_Expected_realized_VMG[df_Expected_realized_VMG['Rang']<32]['Difference Realized New Expected'], label='ecart à new expected')
plt.legend()
plt.show()

En tracant le graphe des écarts avec les deux méthode la différence est frappante. En effet, là où le premier modele n'était fitter en réalité qu'en moyenne et présentait des écarts importants pour les valeurs excentré, ce nouveau modèle présente des écarts en valeurs plus constant, même si présentant un biais.

In [31]:
df_Expected_realized_VMG.sort_values(by='Difference Realized New Expected', ascending=False)[['Voile', 'Rang', 'VMG realized','Difference Realized New Expected', 'Foil','Année lancement']].reset_index(drop=True)
Out[31]:
Voile Rang VMG realized Difference Realized New Expected Foil Année lancement
0 59 2 12.136220 4.091760 True 2019
1 1 3 11.977165 3.428299 False 2007
2 27 14 11.462992 3.362313 True 2007
3 53 13 11.455118 3.226558 False 2007
4 18 6 12.208661 3.095952 True 2015
5 34 10 11.749606 3.095230 True 2015
6 1000 9 11.843307 3.089525 False 2008
7 30 17 11.088189 2.831420 False 2011
8 49 16 11.195276 2.810963 False 2007
9 9 12 11.653543 2.769417 False 2007
10 7 18 10.711024 2.548276 True 2007
11 109 11 11.697638 2.412017 True 2010
12 17 5 12.122047 2.188401 True 2015
13 69 30 7.408661 2.155081 False 2005
14 83 26 7.825984 2.074978 False 2006
15 777 24 8.382677 1.969646 False 1999
16 79 1 12.892063 1.912980 True 2019
17 10 8 12.145669 1.892578 True 2015
18 4 7 12.233071 1.785951 True 2019
19 72 27 7.803150 1.541972 False 1998
20 99 15 11.381890 1.482222 True 2019
21 85 4 12.152756 1.450128 True 2010
22 50 28 7.694488 1.415646 False 2006
23 2 23 8.653543 1.222875 True 2020
24 222 29 7.403150 1.125994 False 2007
25 33 22 8.570079 0.779836 False 2000
26 14 20 8.732258 0.414246 True 2007
27 71 21 8.639370 0.274681 False 2007
28 56 25 8.077165 0.218877 True 2015
29 11 31 7.084252 0.076459 True 2019
30 92 19 9.934646 -0.006813 True 2007
31 8 32 6.189764 -2.469821 True 2018

Ce tableau nous permet ainsi d'obtenir un nouveau classement, qui permettrait d'être plus juste en fonction des bateaux, et donc de comparer les compétences des skippers.

Pour conclure, ce genre de méthode permet en effet de bien rendre compte de l'impact des bateaux dans le classement de la course, ce qui n'est pas très surprennant. Il pourrait ici paraitre attrayant de corriger comme il a été fait dans ce projet les classements en fonction des performances des bateaux. Toutefois il se pose plusieurs limites à cela. D'une part la fiabilité de ces modèles, en effet dans le cas présent étant donné du manque de données (seulement 90 points par bateaux actuellement) le modèle représentant la vitesse de ceux-ci en fonction du vent reste peu fiable et très imparfait. De plus parmi les compétences des skippers le choix et/ou la création de leur navire est une composante à prendre en compte. En effet ils sont, avec leurs sponsors évidemment, les décideurs de leur navire et s'impliquent souvent dans le choix, la création et les modifications de ceux-ci.

C - Utilisation des modèles créés précédement et des données de vent pour prédire le meilleur trajet pour chaque bateau

Ayant désormais des modèles, bien qu'imparfaits, permettant d'estimer la vitesse des bateaux selon l'angle et la vitesse du vent, il est désormais possible d'estimer la route la plus rapide pour finir la course. Pour ce faire j'ai créé les fonctions ci-dessous. Le fonctionnement est assez simple, ayant mis des points balises sur la carte, une fonction va chercher le chemin le plus rapide d'un point à l'autre. Pour ce faire deux méthodes sont utilisées ici. La premiere est la fonction from_a_to_b qui va chercher à estimer le chemin le plus rapide en regardant toutes les 6h quel cap permet de se rapprocher le plus du point objectif en fonction des conditions de vent, en essayant un nombre de déviations défini. La seconde from_a_to_b_with_random_variation, va découper le parcours entre les balises en sous parcours en ajoutant des points incluant des déviations plus ou moins importantes. Cette méthode est appliquée de facon récursive et fait au final appel à la méthode from_a_to_b puis compare les longueurs totale des chemins pour finalement decider quel chemin est le plus rapide. Pour faire tout cela il faut évidement disposer de données de vent, hors les valeurs futures ne sont bien évidemment pas disponibles. Pour gérer cela, bien que de façon imparfaite nous utiliserons les données de vent à la même date et heure mais à l'année passée.

In [32]:
global cap_bonne_esperance_pos
global Leuuwin_pos
global canary
global portugal
global horn
global goal

cap_bonne_esperance_pos=(-39,17.5)
Leuuwin_pos=(-45,128)
horn=(-58,-65)
canary=(20,-26)
portugal=(44,-10)
goal=(46.26,-2.5)
In [33]:
def find_best_way_to_goal(cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,start_date,boats_model, boatnumber,
                          nb_random_mid_point, start_post, next_point=0,nb_dev=1):
    #main function to search the path from the boat to the goal
    boat_speed_model=boats_model.get(boatnumber)
    objectif_list=[cap_bonne_esperance_pos,Leuuwin_pos,horn,canary,portugal,goal]
    bestroad=[]
    if nb_dev>5:
        print("maximum 5 deviation")
        return None
    else:
        for goal_number in range(next_point, len(objectif_list)):
            if goal_number>3:
                spread_size=0.6
            else: 
                spread_size=4
            road,reach_goal=from_a_to_b_with_random_variation(start_post,objectif_list[goal_number],cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,start_date,boat_speed_model,nb_random_mid_point,spread_size=spread_size, nb_dev=nb_dev, first=True)
            if not reach_goal:
                road,_=from_a_to_b_with_random_variation(start_post,objectif_list[goal_number],cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,start_date,boat_speed_model,nb_random_mid_point,spread_size=spread_size, nb_dev=nb_dev, first=True,use_deg_special_untrapper=True)
            bestroad+=road
            start_post=bestroad[-1]
            #print(start_post)
        return bestroad
    
    
def from_a_to_b_with_random_variation(a, b, Unc_file_cy, Vnc_file_cy, Unc_file_ly, Vnc_file_ly, start_date,
                                      boat_speed_model, nb_random,spread_size=1, nb_dev=1, first=False, max_len=None,
                                      use_deg_special_untrapper=False):
    #function use to search the path between two balise
    bestroad, not_trapped=from_a_to_b(a,b,cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,start_date,boat_speed_model,spred_want=spread_size, nb_dev=nb_dev,use_deg_special_untrapper=use_deg_special_untrapper)
    changed=False
    reach_goal=True
    if nb_random>0 and len(bestroad)-1>6:

        random_place_on_path=int(np.round(len(bestroad)/2,0))
        
        randoms_mid_points=get_ecart_to_line(a,b,nb_dev)
        for new_point in randoms_mid_points:
            #print(new_point)
            bestroad_to_mid, _=from_a_to_b_with_random_variation(a,new_point,cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,start_date,
                                                              boat_speed_model, nb_random-1,spread_size=spread_size, nb_dev=nb_dev, first=False,
                                                              max_len=len(bestroad),use_deg_special_untrapper=use_deg_special_untrapper)                   
            bestroad_from_mid, reach_goal_from=from_a_to_b_with_random_variation(bestroad_to_mid[-1],b,cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,
                                                                start_date+len(bestroad_to_mid)*6,boat_speed_model, nb_random-1,
                                                                spread_size=spread_size, nb_dev=nb_dev,first=False, max_len=len(bestroad),use_deg_special_untrapper=use_deg_special_untrapper)                            
            new_road=bestroad_to_mid+bestroad_from_mid
            if (len(new_road)<=len(bestroad) and reach_goal_from) or ((not not_trapped) and reach_goal_from):
                bestroad=new_road
                changed=True
                
    if (not not_trapped) and (not changed):
        br=[]
        for item in bestroad:
            if item not in br:
                br.append(item)
        bestroad=copy.deepcopy(br)
        reach_goal=False
    return bestroad, reach_goal   

def transform_a_b_to_same_base(a,b):
    #This function is used when going from a positive longitude to a negative one
    new_a=copy.deepcopy(a)
    new_b=copy.deepcopy(b)
    if abs(a[1]- b[1])>180:
        if a[1]<0:
            new_a=(new_a[0],360+a[1])
        else:
            new_b=(new_b[0],360+b[1])
    return new_a, new_b
                                      
def from_a_to_b(a, b, Unc_file_cy, Vnc_file_cy, Unc_file_ly, Vnc_file_ly, start_date,boat_speed_model, spred_want=6, nb_dev=1, max_len=None,use_deg_special_untrapper=False):
    #search the best way to go from a to b using only the current wind data to search which option is the best at each date
    path=[]
    path.append(a)
    
    #get lon, lat, and time in file. Reindex the ones in the last year to make them be used while don't have anymore data
    lon_in_file,lat_in_file, time_in_file_cy=get_lon_lat_time_from_nc(cy_Vwind)
    time_in_file_cy-=time_in_file_cy[0]
    _, _, time_in_file_ly=get_lon_lat_time_from_nc(ly_Uwind)
    time_in_file_ly-=time_in_file_ly[0]
    for t in range(len(time_in_file_cy)):
        time_in_file_ly[t]+=max(time_in_file_cy)
    
    if max_len==None:
        max_len=1000
    quarter_day=0
    not_trapped=True
    while get_distance(a,b)>spred_want and quarter_day<max_len and not_trapped:

        windspeed, windangle=get_values_wind(Unc_file_ly, Vnc_file_ly, lat_in_file, lon_in_file, time_in_file_ly, a[0], a[1], start_date )

        options=get_options(a,b,windangle, windspeed, boat_speed_model, nb_dev=nb_dev,use_deg_special_untrapper=use_deg_special_untrapper)
        dist_tmp=[]
        new_points=[]
        
        for opt in options:
            endpoint=where_i_arrive(a, opt, b=b)
            if endpoint[0]>-90 and endpoint[0]<90:
                if check_if_line_on_ocean(a, endpoint):
                    new_points.append(endpoint)
                    esb, bsb=transform_a_b_to_same_base(endpoint, b)
                    dist_tmp.append(get_distance(esb, bsb))
        if len(new_points)>0:
            a=new_points[np.argmin(dist_tmp)]
            path.append(a)
            start_date+=6
            quarter_day+=1
        else: 
            not_trapped=False
            path+=[a for iter in range(1000)]
    asb,bsb=transform_a_b_to_same_base(a,b)
    if  get_distance(asb,bsb)>spred_want:
        not_trapped=False
    return path, not_trapped
    
def get_options(a,b,windangle, windspeed, boat_speed_model, nb_dev=1,use_deg_special_untrapper=False):
    #get the options the boat have at each date
    a, b = transform_a_b_to_same_base(a, b)
    base_vect=np.array((b[0]-a[0],b[1]-a[1]))/get_distance(a, b)
    val=np.arccos(base_vect[0])
    dev=np.pi/2/(nb_dev+1)
    if use_deg_special_untrapper:
        dev=np.pi/(nb_dev+1)
    options=[]
    for dev_i in np.arange(-nb_dev, nb_dev+1,1):
        deg=get_angle_wind_boat(windangle,math.degrees(val+dev*dev_i))
        #print("speed: "+str(windspeed)+" and deg: "+str(deg))
        boat_speed=boat_speed_model.predict(poly.fit_transform([transform_data_for_model(windspeed, deg)[0]])).reshape(-1,1)
        boat_distance_done=boat_speed[0][0]/10
        options.append((np.cos(val+dev*dev_i)*boat_distance_done,np.sin(val+dev*dev_i)*boat_distance_done))
    return options

def check_if_line_on_ocean(a, b):
    is_on_ocean=True
    for iter in range(0,7):
        c=((a[0]*iter+b[0]*(6-iter))/6,(a[1]*iter+b[1]*(6-iter))/6)
        if not globe.is_ocean(c[0], c[1]):
            is_on_ocean=False
    return is_on_ocean

def get_ecart_to_line(a_r, b_r, nb_ecart):
    #create deviation from the most straight road to force the function to test other way of reaching point
    if np.sign(b_r[1])!=np.sign(a_r[1]):
        if a_r[1]>0:
            a,b= a_r, (b_r[0],360+b_r[1])
        else:
            a,b=(a_r[0],360+a_r[1]), b_r
        transformed=True
    else: 
        a,b=a_r, b_r
        transformed=False
    points=[]
    vect_a_b=(b[0]-a[0],b[1]-a[1])
    D=(a[0]+vect_a_b[0]/2,a[1]+vect_a_b[1]/2)
    AD_dist=get_distance(a,D)
    deg=np.sin(np.pi/4)
    C=((a[0]+b[0])/2,(a[1]+b[1])/2)
    if transformed and C[1]>180:
        C=(C[0],C[1]-360)
    points.append(C)
    for deg_num in range(nb_ecart):
        deg=np.pi/3/7*(deg_num+1)
        DC_dist=AD_dist*np.sin(deg)/np.sqrt(1-np.sin(deg)**2)
        DC_vect=(vect_a_b[1]*DC_dist/AD_dist,vect_a_b[0]/2*DC_dist/AD_dist) 
        C=(D[0]+DC_vect[0], D[1]+DC_vect[1])
        if transformed and C[1]>180:
            C=(C[0],C[1]-360)
        if C[0]>-90 and C[0]<90 and C[1]>-180 and C[1]<-180:
            if globe.is_ocean(C[0],C[1]):
                points.append(C)
            else:
                deg=np.pi/3/7/14*(deg_num+1)
                DC_dist=AD_dist*np.sin(deg)/np.sqrt(1-np.sin(deg)**2)
                DC_vect=(vect_a_b[1]*DC_dist/AD_dist,vect_a_b[0]/2*DC_dist/AD_dist) 
                C=(D[0]+DC_vect[0], D[1]+DC_vect[1])
                if transformed and C[1]>180:
                    C=(C[0],C[1]-360)
                if C[0]>-90 and C[0]<90 and C[1]>-180 and C[1]<-180:
                    if globe.is_ocean(C[0],C[1]):
                        points.append(C)
        deg=-np.pi/3/7*(deg_num+1)
        DC_dist=AD_dist*np.sin(deg)/np.sqrt(1-np.sin(deg)**2)
        DC_vect=(vect_a_b[1]*DC_dist/AD_dist,vect_a_b[0]/2*DC_dist/AD_dist) 
        C=(D[0]+DC_vect[0], D[1]+DC_vect[1])
        if transformed and C[1]>180:
            C=(C[0],C[1]-360)
        if C[0]>-90 and C[0]<90 and C[1]>-180 and C[1]<-180:
            if globe.is_ocean(C[0],C[1]):
                points.append(C)
            else:
                deg=-np.pi/3/14*(deg_num+1)
                DC_dist=AD_dist*np.sin(deg)/np.sqrt(1-np.sin(deg)**2)
                DC_vect=(vect_a_b[1]*DC_dist/AD_dist,vect_a_b[0]/2*DC_dist/AD_dist) 
                C=(D[0]+DC_vect[0], D[1]+DC_vect[1])
                if transformed and C[1]>180:
                    C=(C[0],C[1]-360)
                if C[0]>-90 and C[0]<90 and C[1]>-180 and C[1]<-180:
                    if globe.is_ocean(C[0],C[1]):
                        points.append(C)
    return points
    
def where_i_arriveold(start, move, b=None):
    #function to calculate the reach point from a start point and a move calculated before with wind and boat datas
    if b==None:
        new_lon=start[1]+move[1]
        if new_lon>=180:
            new_lon-=360
        elif new_lon<=-180:
            new_lon+=360
        reach_point=tuple([start[0]+move[0], new_lon])
    else:
        if b[1]-start[1]<move[1] and b[0]-start[0]<move[0]:
            reach_point=b
        else:
            new_lon=start[1]+move[1]
            if new_lon>=180:
                new_lon-=360
            elif new_lon<=-180:
                new_lon+=360
            reach_point=tuple([start[0]+move[0], new_lon])

    return reach_point

def where_i_arrive(start, move, b=None):
    #function to calculate the reach point from a start point and a move calculated before with wind and boat datas
    
    new_lon=start[1]+move[1]
    if new_lon>=180:
        new_lon-=360
    elif new_lon<=-180:
        new_lon+=360
    reach_point=tuple([start[0]+move[0], new_lon])
    if b!=None:
        esb, bsb=transform_a_b_to_same_base(start, b)
        esn, bsn=transform_a_b_to_same_base(start, reach_point)
        if get_distance(esb, bsb)<get_distance(esn, bsn):
            reach_point=b
    return reach_point

def get_distance(a, b):
    #return the euclidian distance between two points
    return np.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2)
    


def get_angle_wind_boat(windangle, degree):
    angle=abs(windangle-degree)
    if angle>180:
        angle-=90
    if angle>180:
        angle=angle-180
    return(abs(angle-180))
    

Désormais nous pouvons pour chaque bateau calculé un chemin théorique pour atteindre l'arrivé en utilisant ces fonctions

In [34]:
####Choose parameters here ### Caution, high value mean very long calculation time ###
nb_dev=5
nb_ecart=5

#############################################################################
def modify_values_chemin_for_graph(chemin):
    pass_max_lon=False
    for pos in range(1,len(chemin)):
        if not pass_max_lon and chemin[pos-1][1]>0 and chemin[pos][1]<0:
            pass_max_lon=True
        if pass_max_lon:
            chemin[pos]=(chemin[pos][0],chemin[pos][1]+360)
    return chemin[1:]

if launch_long_calculations:
    df_with_new_values=copy.deepcopy(df_complete)
    df_with_new_values=df_with_new_values[['Voile', 'Skipper', 'position','date']]
    df_with_new_values['date']=df_with_new_values['date'].map(transform_date)
    start_date=np.sort(df_with_new_values['date'].unique())[-1]
    for boat in boats_numbers:
        tmp=df_with_new_values[df_with_new_values['date']==start_date]
        start_pos=tmp[tmp['Voile']==boat]['position'].values[0]
        skipper=tmp[tmp['Voile']==boat]['Skipper'].values[0]
        chemin=find_best_way_to_goal(cy_Uwind,cy_Vwind,ly_Uwind,ly_Vwind,start_date,boats_model_lasso, boat,
                              nb_ecart, start_pos, next_point=0,nb_dev=3)
        chemin=modify_values_chemin_for_graph(chemin)
        chemin_date=[start_date+6*iter for iter in range(1,len(chemin)+1)]
        Voile_current_chemin=[boat for iter in range(len(chemin))]
        Skipper=[skipper for iter in range(len(chemin))]
        df_tmp=pd.DataFrame({'Voile': Voile_current_chemin, 'Skipper': skipper, 'position': chemin, 'date': chemin_date})
        df_with_new_values=df_with_new_values.append(df_tmp)

    delete_if_exist_and_recreate_directory(os.getcwd()+"\\temporary_save_file\\")
    df_with_new_values.to_csv(os.getcwd()+"\\temporary_save_file\\result_to_end_run.csv")
else:
    df_with_new_values=pd.read_csv(os.getcwd()+"\\temporary_save_file\\result_to_end_run.csv")
    df_with_new_values['position']=df_with_new_values['position'].map(lambda x: eval(x))

Et sortir un classement final

In [35]:
tmp=df_with_new_values.groupby(['Skipper']).max().sort_values(by='date')
tmp=tmp.join(df_complete.groupby(['Skipper'])['Foil'].unique().map(lambda x: x[0]))
tmp=tmp.join(df_complete.groupby(['Skipper'])['Année lancement'].unique().map(lambda x: x[0]))
tmp=tmp.join(df_complete[df_complete['date']==df_complete.groupby(['Skipper']).max()['date'].values[0]].groupby(['Skipper']).max()['Rang'])
tmp=tmp.rename(columns={'Rang': 'Dernier Classement'})
tmp['durée parcours (jours)']=tmp['date'].map(lambda x: (x-7502)/24)
tmp=tmp.reset_index().reset_index().rename(columns={"index": "Rang simulé"})[['Skipper','Voile','Rang simulé','Dernier Classement','Foil','Année lancement','durée parcours (jours)']]
tmp['Rang simulé']=tmp['Rang simulé'].map(lambda x: x+1)
tmp
Out[35]:
Skipper Voile Rang simulé Dernier Classement Foil Année lancement durée parcours (jours)
0 Charlie Dalin 79 1 1 True 2019 59.625
1 Sébastien Simon 4 2 7 True 2019 60.125
2 Kevin Escoffier 85 3 4 True 2010 60.125
3 Boris Herrmann 10 4 8 True 2015 61.125
4 Alex Thomson 99 5 15 True 2019 62.125
5 Yannick Bestaven 17 6 5 True 2015 62.375
6 Samantha Davies 109 7 11 True 2010 64.625
7 Louis Burton 18 8 6 True 2015 64.875
8 Damien Seguin 1000 9 9 False 2008 66.125
9 Benjamin Dutreux 9 10 12 False 2007 66.375
10 Jean Le Cam 1 11 3 False 2007 66.875
11 Giancarlo Pedote 34 12 10 True 2015 67.375
12 Isabelle Joschke 27 13 14 True 2007 68.125
13 Stéphane Le Diraison 92 14 19 True 2007 68.875
14 Romain Attanasio 49 15 16 False 2007 69.125
15 Thomas Ruyant 59 16 2 True 2019 69.375
16 Maxime Sorel 53 17 13 False 2007 70.125
17 Clarisse Cremer 30 18 17 False 2011 70.125
18 Alan Roura 7 19 18 True 2007 70.875
19 Arnaud Boissieres 14 20 20 True 2007 76.875
20 Manuel Cousin 71 21 21 False 2007 77.375
21 Armel Tripon 2 22 23 True 2020 79.625
22 Didac Costa 33 23 22 False 2000 79.625
23 Jérémie Beyou 8 24 32 True 2018 79.625
24 Fabrice Amedeo 56 25 25 True 2015 80.125
25 Kojiro Shiraishi 11 26 31 True 2019 84.625
26 Ari Huusela 222 27 29 False 2007 85.375
27 Alexia Barrier 72 28 27 False 1998 86.125
28 Miranda Merron 50 29 28 False 2006 86.375
29 Pip Hare 777 30 24 False 1999 86.625
30 Clément Giraud 83 31 26 False 2006 94.375
31 Sébastien Destremau 69 32 30 False 2005 101.375

On voit que le modèle est quelque peu optimiste avec près de 20 concurent dépassant le record du Vendée Globe ! Cela est explicable par le fait qu'il ne prend pas en compte les autres conditions météorologique, et que la vitesse des bateaux est croissante de la vitesse des vents, il cherchera donc a trouver les zones de vents les plus rapide. Dans la réalité les skippers feront parfois le choix de ne pas prendre les zones les plus venteuse du fait des dangers que cela peut représenter pour le navire et eux même. Enfin il ne prend pas en compte les avaris possible, mais surtout utilise des données de vent fixé à l'avance. Dans la réalité les prévisions s'étalent au maximum sur les 16 prochains jours et il est donc difficiles de comparer une route tracer dans ces conditions à celle de notre algorithme. Il sera toutefois interressant de regarder les trajets emprunté dans la partie suivante.

Partie Visualisation:

Afin de pouvoir visualiser au mieux la course j'ai choisi d'utiliser un graphe animé à l'aide de plotly express. Ce package permet a la fois d'inserer simplement les données sur la carte, mais permet aussi d'obtenir un graphe interactif contenant les informations des bateaux, et d'animer celui-ci.

In [36]:
#Prepare a dataframe for the graph
df_for_graph_tmp=copy.deepcopy(df_complete)
df_for_graph_tmp['lon']=df_for_graph_tmp['position'].map(lambda x: x[1])
df_for_graph_tmp['lat']=df_for_graph_tmp['position'].map(lambda x: x[0])
dates=np.sort(df_for_graph_tmp['date'].unique())
df_for_graph=df_for_graph_tmp[df_for_graph_tmp['date']==dates[0]]
df_for_graph['time_index']=0
for iter in range(1,len(dates)):
    dt=dates[:iter]
    df_tmp=df_for_graph_tmp[df_for_graph_tmp['date'].map(lambda x: True if x in dt else False)]
    df_tmp['time_index']=iter/4
    df_for_graph=df_for_graph.append(df_tmp)

Creation du graphe sur les données réelles

In [37]:
#Create the graph, get patient can take up to one minute because of the animation
    
fig = px.line_mapbox(df_for_graph,
                   lat='lat',
                   lon='lon',
                   color="Skipper",
                   animation_frame = 'time_index',
                   hover_data=['time_index','Voile'],
                   zoom=1,
                   height=800,
                   title = 'Visualization of the Vendee Globe par Skypper')

fig.update_layout(
    mapbox_style="white-bg",
    mapbox_layers=[
        {
            "below": 'traces',
            "sourcetype": "raster",
            "sourceattribution": "United States Geological Survey",
            "source": [
                "https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"
            ]
        }
      ],
    mapbox_center_lat=0,
    margin={"r":0,"t":0,"l":0,"b":0},
)

fig.update_layout()
fig.show()

Creation d'un graphe statique sur les données estimées

Ayant estimé des trajets pour les bateaux dans la partie précédente, en voici la représentation. Ce graphe ci est statique, en effet la quantité de donnée étant beaucoup plus grande, la durée de la création du graphe dynamique est exponentiel avec le nombre d'observation par bateau. Un graphe dynamique est proposé ci-dessous mais attention celui-ci prend un temps très long à se créer.

In [48]:
#Create the graph, get patient can take up to one minute because of the animation
df_for_graph_tmp=copy.deepcopy(df_with_new_values).sort_values(by='date')
df_for_graph_tmp=df_for_graph_tmp.reset_index(drop=True)
df_for_graph_tmp['lon']=df_for_graph_tmp['position'].map(lambda x: x[1])
#df_for_graph_tmp['lon'][df_for_graph_tmp['date']>8258]=df_for_graph_tmp['lon'][df_for_graph_tmp['date']>8258].map(lambda x: 360+x)
df_for_graph_tmp['lat']=df_for_graph_tmp['position'].map(lambda x: x[0])

fig = px.line_mapbox(df_for_graph_tmp,
                   lat='lat',
                   lon='lon',
                   color="Skipper",
                   hover_data=['date','Voile'],
                   zoom=1,
                   height=800,
                   title = 'Visualization of the Vendee Globe par Skypper')

fig.update_layout(
    mapbox_style="white-bg",
    mapbox_layers=[
        {
            "below": 'traces',
            "sourcetype": "raster",
            "sourceattribution": "United States Geological Survey",
            "source": [
                "https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"
            ]
        }
      ],
    mapbox_center_lat=0,
    margin={"r":0,"t":0,"l":0,"b":0},
)

fig.update_layout()
fig.show()

On peut observer que tout les bateaux n'empruntent pas le même chemin et ce du fait des positions de départ différentes, mais aussi des modèles de vitesse prédite en fonction du vent propre à chacun. Au final l'algorithme de routage fonctionne bien, l'ensemble des navires arrivant à trouver un chemin jusqu'à destination et ce sans se coincer dans un obstacle (ce qui n'a pas toujours été aisé à configurer de prime abord).

Creation d'un graphe Dynamique sur les données estimées

Cette partie n'est lancé que si vous choisissez launch long calculation,le temps estimé de création est de 1h. Attention !!! Si lancé, le résultat sera bien visible mais le fichier ne sera plus enregistrable. Si celui-ci est enregistré avec le graphe alors il ne sera plus possible d'ouvrir le fichier car il causera une erreur mémoire.

In [ ]:
#Prepare a dataframe for the graph
if launch_long_calculations:
    df_for_graph_tmp=copy.deepcopy(df_with_new_values).sort_values(by='date')
    df_for_graph_tmp=df_for_graph_tmp.reset_index(drop=True)
    df_for_graph_tmp['lon']=df_for_graph_tmp['position'].map(lambda x: x[1])
    df_for_graph_tmp['lat']=df_for_graph_tmp['position'].map(lambda x: x[0])
    dates=np.sort(df_for_graph_tmp['date'].unique())
    df_for_graph=df_for_graph_tmp[df_for_graph_tmp['date']==dates[0]]
    df_for_graph['time_index']=0
    for iter in range(1,len(dates)):
        dt=dates[:iter]
        df_tmp=df_for_graph_tmp[df_for_graph_tmp['date'].map(lambda x: True if x in dt else False)]
        df_tmp['time_index']=iter/4
        df_for_graph=df_for_graph.append(df_tmp)
    df_for_graph.to_csv(os.getcwd()+"\\temporary_save_file\\result_to_end_run_for_graph.csv")
else:
    df_for_graph=pd.read_csv(os.getcwd()+"\\temporary_save_file\\result_to_end_run_for_graph.csv")
In [ ]:
#Create the graph, get patient can take up to one hour because of the animation
if launch_last_graph:   
    fig = px.line_mapbox(df_for_graph,
                       lat='lat',
                       lon='lon',
                       color="Skipper",
                       animation_frame = 'time_index',
                       hover_data=['time_index','Voile'],
                       zoom=1,
                       height=800,
                       title = 'Visualization of the Vendee Globe par Skypper')

    fig.update_layout(
        mapbox_style="white-bg",
        mapbox_layers=[
            {
                "below": 'traces',
                "sourcetype": "raster",
                "sourceattribution": "United States Geological Survey",
                "source": [
                    "https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"
                ]
            }
          ],
        mapbox_center_lat=0,
        margin={"r":0,"t":0,"l":0,"b":0},
    )

    fig.update_layout()
    fig.show()